

clear


%% data
load('data_5y.mat')
d1=data(isfinite(data.y)&isfinite(data.D),:);
d1(d1.year>=1914&d1.year<=1918,:)=[];
d1(d1.year>=1939&d1.year<=1945,:)=[];
ccodes=unique(d1.ccode);
years=unique(d1.year);
N=length(ccodes);
T=length(years);
y=zeros(N,T); D=y; M=y;
for n=1:N
    dn=d1(d1.ccode==ccodes(n),:);
    for t=1:T
        if sum(years(t)==dn.year)>0
            if isfinite(dn.y(years(t)==dn.year))&&isfinite(dn.D(years(t)==dn.year))
                M(n,t)=1; % not missing
                y(n,t)=dn.y(years(t)==dn.year);
                D(n,t)=dn.D(years(t)==dn.year);
            end
        end
    end
end
clear d1 dn n t

dd=readtable('../../../Data/distance_capitals.csv');
Z=10^(-6)*table2array(dd(:,2:end));


%% for Knitro
options=optimset('Display','off');


%% starting values

% "short" replication
load('tDGProbustTH_replication.mat','DGP')
pars0=DGP;

% % "long" replication
% load('tDGProbustTH_replication.mat','srng')
% rng(srng)
% feas=0; % ensure finite log_posterior starting value
% while feas==0
%     beta=[randn(N,1);randn(N,1)]/100;
%     s=rand(N,1);
%     ab=2*rand;
%     lambda=2*rand;
%     a=2*rand;
%     xi=2*rand;
%     b0A=randn/100;
%     b0D=randn/100;
%     pars0=[beta;s;ab;lambda;a;xi;b0A;b0D];
%     feas=isfinite(log_posterior_trueDGP(pars0,y,D,M,Z));
% end 


%% parameter estimates
[DGP,lp,exit]=knitromatlab(@(x)log_posterior_trueDGP(x,y,D,M,Z),pars0,[],[],...
  [],[],[-Inf(2*N,1);zeros(N+4,1);-Inf(2,1)],Inf(length(pars0),1),[],[],...
  options,'knitro.opt');


%% growth shocks variance matrix
load('prior_robustTH.mat')

s=DGP(2*N+(1:N));
a=DGP(3*N+3);
xi=DGP(3*N+4);

Sigma=diag(s)*(a*exp(-xi*Z))*diag(s);
Sigma=tril(Sigma)+tril(Sigma,-1)';  % enforce exact symmetry

clearvars -except a0 b0A b0D d_v d_vs f0 om_a om_b om_f om_th ...
    s_v s_vs th0 Sigma

save('priorWSigma_robustTH.mat')



